The importance of the electronic contribution to linear magnetoelectricity 
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We demonstrate that the electronic contribution to the linear magnetoelectric response, usually 
omitted in first-principles studies, can be comparable in magnitude to that mediated by lattice 
distortions, even for materials in which responses are strong. Using a self-consistent Zeeman response 
to an applied magnetic field for noncollinear electron spins, we show how electric polarization emerges 
in linear magnetoelectrics through both electronic- and lattice-mediated components - in analogy 
with the high- and low-frequency dielectric response to an electric field. The approach we use is 
conceptually and computationally simple, and can be applied to study both linear and non-linear 
responses to magnetic fields. 

PACS numbers: 75.85.+t, 71.15.Mb, 75.30.Cr, 



Linear magnetoelectrics are materials which respond 
with a change in electric polarization to a magnetic field, 
or conversely with a change in magnetization to an elec- 
tric field: Pi = ouaHa'. Ma 



ctijEi 



where a is the 



linear magnetoelectric tensor. 

The research challenges in identifying materials with 
useful linear magnetoelectric (ME) responses are three- 
fold: (i) Symmetry requirements that both space- 
inversion and time-reversal symmetries be broken are sat- 
isfied by few materials, (ii) materials satisfying these cri- 
teria tend to do so only in phases that develop at rela- 
tively low temperatures, and (iii) most of the MEs dis- 
covered to date have weak responses. Recently, a num- 
ber of developments have led to a significant revival of 
activity in the search for novel magnetoelectric materi- 
als, including the observation that multiferroics can have 
strong ME responses [l], and the identification of the mi- 
croscopic couplings that are responsible for strong and 
weak ME responses 

First-principles methods are emerging d, 0] as a valu- 
able tool for computing the strength of a in real ma- 
terials without any empirical input. The methods are 
becoming sufficiently reliable to be used in a predictive 
capacity in searching for new ME materials. While many 
approaches have been explored or can be envisaged for 
computing a, ranging from a self-consistently applied 
electric or magnetic field, to quantum-mechanical pertur- 
bation theory in the applied field, by far the most suc- 
cessful and widely used approach for applications to date 
has been a linear-response approach based on the lattice- 
dynamical quantities [6|-|8|. However, this approach com- 
putes only the so-called "lattice-mediated" part of a and 
ignores purely electronic contributions to the response. 
The common justification is that such contributions are 
expected to be weak, just as in strong dielectrics the elec- 
tronic response is negligible compared to the ionic con- 
tribution. In this Letter, we demonstrate, using an alter- 
native numerical approach involving a self-consistently 
applied magnetic field, that the purely electronic magne- 



toelectric response can in fact be large, even in materials 
with relatively strong a in which one might expect the 
response to be dominated by lattice mechanisms. In fact, 
the electronic and ionic contributions can be of similar 
magnitude and opposite sign, which can lead to a weak 
total response while lattice-only methods would predict 
it to be large. 

We begin by defining the different contributions to a. 
We define a "clamped-ion" contribution that accounts for 
magnetoelectric effects occuring in an applied field with 
all ionic degrees of freedom remaining frozen. This purely 
electronic contribution, a cl , is the response that would be 
measured for high-frequency fields, in analogy with the 
static-high-frequency dielectric response in insulators, 
too. The remaining part of the response, which emerges 
at low frequency so that the ions (and in principle the 
lattice parameters respond to the field, we label a latt . 
a iatt j g |- ne difference between the total response and the 
electronic, so that a tot = a cl + a latt , and it is the part 
of the response that is accessible from linear-response 
theory using lattice-dynamical quantities (Note that 
Q iatt a j so j nc i U( ies some electronic response through the 
dynamical charges). 

Evaluating a o1 requires a full quantum-mechanical 
treatment of the response, either using perturbation the- 
ory or a self-consistently-applied finite field. In this work, 
we use an applied magnetic field to compute a. Due to 
difficulties with periodic boundary conditions when us- 
ing a full vector potential in the electronic Hamiltonian, 
we restrict our field to act only on the electron spin as 
a self-consistent Zeeman field. This means that we omit 
contributions to a that are derived from the orbital mag- 
netic response[10], which are expected to be significantly 
weaker than the spin-derived response for most systems. 

Computational Approach: Calculations of ME re- 
sponses require a treatment of noncollinear spin or- 
ders. Within the Kohn-Sham framework, noncollinear- 
ity is handled by generalizing the orbitals to be complex 
spinors, resulting in a 2 x 2 spin-density matrix (n CT(T ') 
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that allows the magnetization density to vary in both 
magnitude and direction throughout the system. 

In order to apply a magnetic field, we begin by 
making a Legendre transform of the noncollinear-spin 
Hohcnberg-Kohn energy functional: 



A 



(r);H = E HK [n aa > (r)] - fi H ./2 to 



(1) 



where £hk is the usual zero-field energy functional (con- 
sisting of non-interacting kinetic energy, external elec- 
trostatic energy, Hartree and exchange-correlation terms, 
and the Madelung energy of ion-ion interactions). H is 
the auxiliary magnetic field applied to spin degrees of 
freedom and /2 to t is the total spin magnetic moment of 
the system (including the Lande g factor). The spatially 
varying magnetic moment, fi(r), can be found from the 
four components of n aa > , and /2 to t is its spatial integral. 
Subsequently, we variationally minimize A with respect 
to single-particle Kohn-Sham spinor orbitals, subject to 
the usual constraints on orthonormality and conserva- 
tion of total particle number. The result is a term in the 
Kohn-Sham potential, acting on the spinor orbitals, that 
imposes the external magnetic field on the noncollinear 
spin density. Practically, this term simply shifts the rela- 
tive external potential for each of the four spin manifolds: 
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(2) 



which is trivially compatible with periodic boundary con- 
ditions, and which clearly reduces to the collinear case, 
with the field providing a different Fermi level for "up" 
and "down" spin channels, if H = (0,0, H z ). 

We implement Eq. p into the Vienna Abinitio Simula- 
tion Package (VASP)[11(. We employ a plane- wave basis 
set for expanding the electronic wave functions and den- 
sity, and PAW potentials [l2| are used for core- valence 
separation. We note that it is important to disable sym- 
metrization of the wave functions in the Brillouin zone 
since application of a magnetic field leads to an elec- 
tronic structure that breaks the crystal symmetry. We 
self-consistently include spin-orbit coupling in all of our 
calculations. This is required to obtain a magnetoelec- 
tric response for those materials in which the spin-lattice 
coupling derives from relativistic effects such as the an- 
tisymmetric (s*j x Sj ) Dzyaloshinskii-Moriya interaction. 
All of our calculations simulate a single magnetic domain, 
which corresponds to experiments in which a poling pro- 
cedure has been performed. Since different antiferromag- 
netic domains contribute to a with different signs, mea- 
surements otherwise tend to provide a lower bound on 
the single-domain response. 

Transverse magnetoelectric response of Cr%0^: 
Cr2C>3, the first ME material to be discovered, remains 
the best-studied and prototypical linear magnetoelectric. 
Cr 2 C>3 adopts the space group R3c in the ground state, 
with Cr and O occupying Wyckoff positions 4c and 



6e in the rhombohedral setting. We work with the 
experimental volume (95. 9A 3 ) and rhombohedral angle 
(55.13°)[T3], but fully optimize the ion coordinates 
to the LDA+U ground-state (Cr: x = 0.1536; O: 
x = 0.9426). For the exchange-correlation potential, we 
use the Dudarev form of LSDA+J7 with U cB = 2.OeV0j. 
In the ground-state G-type antiferromagnetic (AFM) 
ordering, the free energy of the system contains two 
symmetry allowed linear magnetoelectric couplings and 
is of the formllJI 



F = -a x (E X H X + E y H y ) - a\\E z B z 



(3) 



We note that the response parallel to the trigonal (easy) 
axis, an, has been experimentally demonstrated to be 
close to zero at zero temperature This small low- 
temperature response is expected: In the absence of 
quantum fluctuations and orbital contributions, the low- 
temperature ME response should tend to zero because 
the spin-only parallel magnetic susceptibility vanishes at 
T = OK. Our calculations, which are formally at zero 
temperature and do not include quantum spin fluctua- 
tions or the orbital magnetic response, do give a\\ = 
as expected. The temperature dependence of an has re- 
cently been explained using an alternative first-principles 
scheme [lij. 

For the transverse response of Cr203, a±, we first com- 
pute the electronic ME contribution with clamped ions. 
Upon application of a Zeeman field perpendicular to the 
collinear spin axis, the Cr spins cant so that the unit 
cell acquires a net magnetization. Our calculated spin- 
only magnetic susceptibility is x± = 1-9 x 10~ 3 in di- 
mensionless SI units, which compares favorably to the 
experimental value [r^ of 1.7 x 10~ 3 at 78 K. With the 
new spin configuration generated by the applied magnetic 
field, spin-orbit coupling leads to an electric polarization. 
Note at this stage that the lattice has not been permit- 
ted to relax in response to the field. We compute the 
electric polarization using the Berry phase approach [l8j|. 
The results are shown as open squares in Fig.[TJ with our 
calculated = 0.34ps.m _1 . 

The canting of spins in the applied magnetic field also 
causes spin-orbit-driven forces on ions. The induced 
forces are very small, mirroring the observed weakness 
of the linear ME effect in Cr203, so high numerical qual- 
ity must be achieved in studying the structural distor- 
tion. We carefully converge our ionic forces with re- 
spect to basis-set size (E cut = 700 eV), fe-point sampling 
(6x6x6), and self-consistent field iteration. We then 
choose a small tolerance to which forces are eliminated 
(< 5^eV A" 1 ). 

Once the structural distortions are found for a given 
magnetic-field strength, we quickly estimate the ionic 
polarization by multiplying the ion displacements by 
the full Born-effective-charge tensor (see open circles of 
Fig. [TJ. This procedure leads to the lattice-mediated 
ME response, a 1 ^ = 1.11 ps.m -1 . We also computed 
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Magnetic Field, n H (T) 

FIG. 1. (Color online) The separate contributions to the 
transverse magnetoelectric response of Cr2C>3 calculated us- 
ing Eq. [2] The clamped-ion response, a cl , is shown as open 
squares, and contributes approximately one fourth of the to- 
tal response (filled circles). The remainder of the response is 
due to structural distortions in the applied field (open circles), 
computed using the Born effective charge tensor. 

the same quantity using our parameters with the lattice- 
dynamical method for cv latt introduced previously^ and 
find CJ.9ps.m~ 1 , in reasonably good agreement with the 
applied magnetic field approach. 

With the lattice distortion included for each magnetic- 
field strength, we now fully recompute the electric polar- 
ization using the Berry-phase approach with the magnetic 
field applied. This procedure yields the full spin-mediated 
a, both electronic and lattice contributions, and is shown 
as filled circles in Fig. [TJ We find a± — 1.45 ps.m -1 
(= 4.3 x 10 -4 emu-CGS), in excellent agreement with 
summing a°l_ and erf** as expected. From this analysis, it 
is clear that the electronic response, cv o1 , contributes one 
fourth of the total spin-driven magnetoelectric response 
of Cr2C>3. The total response that we compute is in good 
agreement with zero-temperature extrapolations of ex- 
perimental measurements, which range 0, EH H3] from 
2—4.7 x 1CT 4 emu-CGS. 

Off- diagonal ME response of LiNiPO^: The lithium 
orthophosphates LiMP04 (M=transition metals) are 
currently attracting much interest because of their po- 
tential use as cathodes for Li-ion batteries, as well as 
for their large magnetoelectric responses and the recent 
observation of ferrotoroidic domains in LiCoP04[2l|. 

Here we focus on LiNiP04 as a test case for the size 
of the electronic ME response. LiNiP04 has space group 
Pnma with 28 atoms in the primitive unit cell and four 
Ni 2+ magnetic sites. We again use the experimental [22l| 
lattice parameters for our calculations (a = 10.032 A, 
b = 5.854 A, c = 4.677 A) while the ionic coordinates are 
fully relaxed to the LDA+U ground state. 

We first address the zero-field magnetic structure. 

r 1 

Experimentally 23] , the magnetic structure has been 



characterized as C-type AFM with spins oriented pre- 
dominantly along the c axis, combined with a weak spin 
canting of A-type AFM order along the a axis. The re- 
ported zero-field canting angle, 9, in the Ni compound is 
7.8°. 

With a plane- wave cutoff of 500 eV and 4x4x2 fc-point 
sampling, we qualitatively reproduce the observed mag- 
netic structure in our calculations. However, the precise 
canting angle is strongly sensitive to the intra-atomic J 
parameter used in the Liechtenstein LDA+C/ procedure. 
For U — 5 eV, a common value for Ni 2+ , the canting an- 
gle varies from 1.6° to 7.6° as J is modified from 0.0 
to 1.7 eV. Hence, a relatively large J appears to be im- 
portant for quantitative agreement with the magnitude 
of the reported spin canting. The observed sensitivity 
arises from the double-counting term in the LSDA+U po- 
tential, in which the off-diagonal elements (acting on the 
non-collinear part of the spin density matrix) are deter- 
mined by J alone. The canting angle is quite insensitive 
to our other simulation parameters, including the on-site 
Coulomb interaction, U. 




Magnetic Field, ji H (T) Zero-Field Canting Angle (°) 



FIG. 2. (Color online) L1N1PO4 data under magnetic field H || 
a: (a) Clamped-ion (black) and ionic (red) contributions to 
the polarization versus the magnetic field for different values 
of the intra-atomic Hund's coupling J. (b) Ionic polarization 
versus canting angle for three values of the magnetic field. 



We now turn to the magnetoelectric response of 
LiNiP04. For the magnetic point group mm'm, the mag- 
netoelectric tensor admits only two non-zero components: 
ct xz and ct zx . We focus on the larger component a ZXl 
which we access by computing the Berry-phase polariza- 
tion along c while applying a magnetic field along the a 
axis (H perpendicular to the easy axis). In Fig. [2^l we 
report both the clamped-ion and lattice-mediated mag- 
netoelectric response under this orientation of magnetic 
field. For each of the four different values of the intra- 
atomic J, we find af x ~ l.Ops.m -1 . However, as might 
be expected from the J dependence of the spin canting 
angle, the ionic contributions to a is strongly sensitive to 
J. For J = 0.0 eV, or*^ is moderately large and nega- 
tive. With increasing J, this component of the response 
progressively increases, eventually turning positive. The 
connection between the J dependence of a latt and the 
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FIG. 3. (Color online) Mode contributions to the ionic polar- 
ization of LiNiPCU for a field strength of 20 T. Inset: linear 
evolution of each of the 13 B\ u (polar) lattice modes with 
applied field. 



canting angle is highlighted in Fig. [2b, where the ionic 
contribution to the polarization at different values of the 
magnetic held strength arc plotted against different ini- 
tial canting angles (obtained by varying J). We observe 
a linear relationship between the magnetic-held induced 
ionic polarization and the zero-held canting angle. In- 
terestingly, for a critical canting angle (9 C ~ 3.2°), the 
ionic contribution to the polarization is zero for all held 
strengths. Under these conditions, 0'°* is dominated by 
the electronic contribution, while for other values of the 
canting angle the electronic and ionic components are 
comparable in magnitude. 

To characterize the microscopic origin of the ionic con- 
tribution to a^J* in L1N1PO4, we report in Fig. [3] the 
phonon mode decomposition of the magnetic-field in- 
duced lattice polarization for different values of J. For 
LiNiP04, thirteen polar {B\ u symmetry) modes enter in 
determining a latt . As expected, each mode increases in 
strength linearly with the field (inset of Fig. |3]) , demon- 
strating that the lattice remains harmonic. Plotting the 
con tribution to the polarization for each mode for an ap- 



tions to responses are neglected within this approach, but 
both electronic and ionic contributions to magnetic-field 
induced electric polarization are present. 

We have demonstrated that the electronic contribu- 
tion to magnetoelectric responses is not negligible, and 
can in fact dominate the total response, in contrast with 
prior expectations. Furthermore, we have shown for the 
case of LiNiP04 that the magnetoelectric response is not 
dominated by the softest polar lattice modes, and that 
a partial cancellation of the contributions from different 
modes occurs, weakening the magnetoelectric response of 
this material. These features suggest that the best route 
to engineering strong magnetoelectrics may not lie with 
strain engineering to induce lattice destabilization, but 
rather with strengthening the spin-lattice coupling, for 
example using strong non-relativistic mechanisms [6|. 
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plie A magnetic field of 20 T (Fig. 13")) . wc find that all 



modes are important, so that the response is surpris- 
ingly not dominated by the softest lattice modes. The 
zero a latt at 9 C clearly arises from an accidental cancel- 
lation between the positive and negative magnetoelectric 
responses of the 13 individual polar lattice modes. 

Conclusion: We have demonstrated a new approach 
for calculating the linear magnetoelectric response of ma- 
terials from first principles. The approach, which pro- 
vides an efficient way to extract a even for systems with 
low I symmetry, involves self-consistent application of a 
magnetic field which acts on the spin degrees of freedom 
only in the Zeeman sense. Orbital magnetic contribu- 



[9 
[10 

[11 

[12 
[13 

[14 

M 



116 



Authors contributed equally. 

M. Fiebig, J. Phys. D: Appl. Phys., 38, R123 (2005). 

T. Kimura, T. Goto, H. Shintani, K. fshizaka, T. Arima, 

and Y. Tokura, Nature, 426, 55 (2003). 

T. Goto, T. Kimura, G. Lawes, A. P. Ramirez, and 

Y. Tokura, Phys. Rev. Lett., 92, 257201 (2004). 

H. Katsura, N. Nagaosa, and A. V. Balatsky, Phys. Rev. 

Lett., 95, 057205 (2005). 

S. Picozzi, K. Yamauchi, B. Sanyal, I. A. Sergienko, and 

E. Dagotto, Phys. Rev. Lett., 99, 227201 (2007). 

K. T. Delaney, M. Mostovoy, and N. A. Spaldin, Phys. 

Rev. Lett., 102, 157203 (2009). 

J. Iniguez, Phys. Rev. Lett., 10 1, 117201 (2008). 

J. C. Wojdel and J. Iniguez, Phys. Rev. Lett., 103 



267205 (2009) 



J. C. Wojdel and J. Iniguez, ArXiv:1004.3152vl. 

A. Malashevich, I. Souza, S. Coh, and D. Vanderbilt, 

ArXiv:1002.0300v2. 

G. Kresse and J. Furthmiiller, Phys. Rev. B, 54, 11169 
(1996). 

G. Kresse and D. Joubert, Phys. Rev. B, 59, 1758 (1999). 
L. W. Finger and R. M. Hazen, J. Appl. Phys., 51, 5362 
(1980). 

I. E. Dzyaloshinskii, J. Exptl. Theoret. Phys. (USSR), 
37, 881 (1959). 

V. J. Folen, G. T. Rado, and E. W. Stalder, 



Lett., 6, 607 (1961) . 

M. Mostovoy, A. Scaramucci , 

N. Spaldin, ArXiv:1004.2070vl. 



Phys. Rev 



K. T. Delaney, and 



5 



[17] T. R. McGuire, E. J. Scott, and F. H. Graunis, Phys. 

Rev., 102, 1000 (1956). 
[18] R. D. King-Smith and D. Vanderbilt. Phys. Rev. B. 47, 



1651 (1993). 



[19] H. Wiegelmann, A. G. M. Jansen, P. Wyder, J.-P. Rivera, 

and H. Schmid, Ferroelectrics, 162, 141 (1994). 
[20] E. Kita, K. Siratori, and A. Tasaki, J. Appl. Phys., 50, 
7748 (1979). 



[21] B. B. Van Aken. J.-P. Riv era. H. Schmid, and M. Fiebig, 



Nature, 449, 702 (2007) 



[22] I. Abra hams and K. S. Easson, Acta Cryst. C, 49, 925 



(1993) 



[23] T. Jensen, N. Christensen, M. Kenzelmann, H. R0nnow, 
C. Niedermayer, N. Andersen, K. Lefmann, J. Schefer, 
M. v. Zimmermann, J. Li. J. L. Za restkv, and D. Vaknin, 
Phys. Rev. B, 79, 092412 (2009) 



